#include <it/vec.h>
#include <it/mat.h>
#include <it/source.h>
#include <it/random.h>
#include <it/linalg.h>

#include "utils.h"
#include "constants.h"
#include "upsampling.h"


mat make_carriers( int Nc, int Nv, unsigned long int key, uint ups )
{

  mat temp; 
  mat u; 
  vec v;

  int i= 0; 

  Nv /= ups*ups;
  
  temp = mat_new_zeros( (size_t)Nv, (size_t)Nc );

  it_seed( key );

  for ( i= 0; i< Nc; i++ ) 
    {
      v = source_gaussian( (size_t)Nv, 0.0, 1.0 ); 
      vec_decr( v, vec_mean( v ) ); 
      mat_set_col( temp, (size_t)i, v );
      vec_delete(v);
    }  

  mat_gs( temp ); 

  for ( i= 0; i< Nc; i++ ) 
    {
      v = mat_get_col( temp, i ); 
      vec_div_by( v, sqrt(vec_variance(v)) );

      mat_set_col( temp, i, v );
      vec_delete(v);
    }

  u = mat_upsample_carriers( temp, ups );

  mat_delete( temp );

  return( u );

}

vec get_correlations( mat y, uint key, uint sz, uint Nc, bvec mest, uint ups )
{
  uint Nv = sz*sz;
  mat U = make_carriers( Nc,Nv,key, ups );
  vec corrs = vec_new_zeros( Nc );

  printf("Nc: %d   Nv: %d  sz: %d    mat_height(y): %d      mat_wigth(y): %d \n",Nc,Nv,sz,mat_height(y),mat_width(y));

  int i,j,k,l,nbY=0;
  mat y_fond = mat_new_zeros((size_t)sz,(size_t)sz);
  for(i=0;i<mat_height(y)/sz;i++){
	for(j=0;j<mat_width(y)/sz;j++){	
		nbY++;
  		for(k=0;k<sz;k++){
			for(l=0;l<sz;l++){
			  y_fond[k][l] += y[i*sz+k][j*sz+l];
			}
		}
	}
  }	

  for(k=0;k<mat_height(y_fond);k++)
    for(l=0;l<mat_width( y_fond);l++)
       y_fond[k][l]= y_fond[k][l]/nbY;

  for(i=0;i<Nc;i++)
  {
     vec V = mat_get_col(U,i);
     vec Y = mat_to_vec(y_fond);
     corrs[i] = vec_inner_product(V,Y);
     if(corrs[i]>0) {mest[i]=0;}
	 else {mest[i]=1;}
     vec_delete(Y);
     vec_delete(V);

  }

 mat_delete( U );
 return corrs; 
}

